J1-J2 frustrated two-dimensional Heisenberg model: Random phase approximation 

and functional renormalization group 
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We study the ground state properties of the two-dimensional spin- 1/2 J\- J2-Heisenberg model on 
a square lattice, within diagrammatic approximations using an auxiliary fermion formulation with 
exact projection. In a first approximation we assume a phenomenological width of the pseudofermion 
spectral function to calculate the magnetization, susceptibilities and the spin correlation length 
within RPA, demonstrating the appearance of a paramagnetic phase between the Neel ordered 
and Collinear ordered phases, at sufficiently large pseudo fermion damping. Secondly we use a 
Functional Renormalization Group formulation. We find that the conventional truncation scheme 
omitting three-particle and higher order vertices is not sufficient. We therefore include self-energy 
renormalizations in the single-scale propagator as recently proposed by Katanin, to preserve Ward 
identities in a better way. We find Neel order at g = J2/J1 < g c i ~ 0.4 . . . 0.45 and Collinear order 
at g > g C 2 ~ 0.66 . . . 0.68, which is in good agreement with results obtained by numerical studies. 
In the intervening quantum paramagnetic phase we find enhanced columnar dimer and plaquette 
fluctuations of equal strength. 



I. INTRODUCTION 



It has been known for a long time that quantum an- 
tiferromagnets, i.e. spin-1/2 systems coupled by Heisen- 
berg exchange interaction, are strongly affected by quan- 
tum fluctuations at low temperatures. Thermal fluctu- 
ations are important as well, especially since they sup- 
press long-range-order (LRO) in two dimensions at any 
finite temperature, but their role is relatively well under- 
stood. By contrast, quantum fluctuations operate in a 
much more complex way: they may suppress LRO, but 
may at the same time lead to novel ground states known 
under the labels "spin liquids, valence-bond solids" . The 
first such state proposed in the literature is Anderson's 
RVB-( "resonating- valence-bond" )-statei. In the context 
of cuprate superconductors, viewed as hole-doped Mott 
insulators, RVB-states have been proposed by Anderson 2 - 
to form the fundamental basis on which the theory of 
high-Tc superconductivity should be built. Although the 
idea has been considered by many authors since then, 
there is no conclusive answer to the question of the role 
of a spin liquid state for HTSC. These studies have raised 
the question, however, under which conditions quantum 
fluctuations are strong enough to destroy long range or- 
der. In general, spin liquid type states may be expected 
to be stabilized by any type of quantum fluctuations. 
For spin systems it has been proposed that frustration 
either by competing spin interactions or due to special 
geometric arrangements may lead to a spin liquid state. 
In particular, by tuning the interactions or the lattice 
anisotropy a quantum phase transition from a state with 
long-range order into a spin liquid state may take place. 
Generally speaking it has proven to be, may be unex- 
pectedly, hard to destroy long range order by quantum 
fluctuations. 

The simplest theoretical model of such a system 
is the Quantum Heisenberg Antiferromagnet (QHAF) 
with nearest neighbor interaction for spins h on a two- 



dimensional square lattice. Its ground state is known to 
be the Neel state with staggered magnetization reduced 
by quantum fluctuations^. At any finite temperature the 
magnetic order is destroyed by thermal fluctuations, but 
the correlation length is found to increase exponentially 
with decreasing temperature '~ 5 . This physics has also 
been obtained from the Quantum Nonlinear Sigma Model 
(QNLcrM) in the renormalized classical regime-. 

A simple model with competing interactions is the 
J1-J2 model 7 , featuring an additional antiferromagnetic 
next-nearest neighbor interaction J2 in addition to the 
nearest neighbor coupling J\ . This model has attracted 
attention as a simplified model 8 for the effect of doping 
in the cuprate superconductors: when a small concentra- 
tion of holes is doped into the CuO-planes, the long-range 
AF order of the undoped system is rapidly destroyed*^, 
giving way to a non-magnetic "pseudo-gap" state and to 
superconductivity. 

Recently, this model has also found use for certain 
Vanadate compound a 11 ' 12 , for which the magnetic in- 
teractions can be modeled by the J1-J2 Hamiltonian of 
weakly coupled planes. 

Even more recently the J\-J% model has been invoked 
to account for the weakened AF long range order in the 
iron pnictides^ - — . The universally observed linear tem- 
perature dependence of the magnetic susceptibility of 
these compounds has also been addressed in the frame- 
work of the J1-J2 models. 

If the spins of the model are considered as classical 
(S — > 00), there is an abrupt transition from Neel or- 
der to the Collinear configuration for sufficiently strong 
frustration, J2/J1 = 1/2. In mean field approximation a 
first order transition from the Neel to the Collinear state 
is found. However, this is changed by quantum fluctu- 
ations. Early on it has been found^ ' 17 ' 18 that a non- 
magnetic phase exists in the region 0.4 < J1/J2 < 0.65, 
between the two ordered states. The nature of this in- 
termediate state is what we would like to unravel. Most 
recent mainly numerical work (see Refs. [l9l - |2l1 and refer- 
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ences therein) on the J\-J% model indicates that it may 
be a valence-bond solid 2 ^ (VBS), rather than a homoge- 
neous spin liquid^. In the VBS the spins in the plane 
form pairwise singlets, which are spontaneously dimer- 
ized in a, e.g. columnar pattern and therefore break 
the lattice translational symmetry. It has also been pro- 
posed that the dimerization takes place on units of 2 x 2 
plaquettes 2 ^. Evidence for a VBS has also been found 
in studie d 25 ' 26 of a model of coupled spin chains^ 7 -, when 
the results are extrapolated to the isotropic J\-J% model 
in the plane. Concerning the nature of the quantum 
phase transitions recent studies indicate that the tran- 
sition from the paramagnetic phase to the Collinear con- 
figuration is of first order— ~ 21 ' 28 . On the other hand the 
properties of the transition from the Neel phase to the 
paramagnetic phase are still highly controversial. Recent 
studies point to either a first orde r 19 ' 25 or a second or- 
der transitio n 28 ' 29 . The latter scenario gives rise to the 
question of how two differently ordered phases may be 
connected by a continuous phase transition 3 - .. 

In this paper we develop new analytical methods for 
calculating the ground states and the excitation spectra 
of spin models with competing interactions, such as the 
model discussed above, on the basis of infinite rcsumma- 
tions of perturbation theory in the couplings J\, J^. To 
this end we use a representation of the spin operators 
in terms of pseudofermions^i. One motivation for using 
a fermionic representation rather than a bosonic repre- 
sentation is the available experience in describing spin 
liquids or dimerized spin-singlet states with fermions, 
mainly within large-N and Mean-Field approaches (see 
e.g. Refs. l32l - [35l ). On the other hand, pseudofermion 
representations have hardly been used to study magnetic 
ordering phenomena 36 . Although a large body of results 
of numerical studies of these models is available, ana- 
lytical approaches starting from a microscopic Hamil- 
tonian are rare. We use a newly developed implemen- 
tation of the Functional Renormalization Group (FRG) 
metho d 37 ' 38 applied to interacting quantum spin models. 
In this we are aided by the experience we have previously 
gained with the nearest neighbor Heisenberg mode l 39 ' 40 . 
Auxiliary particle representations of spin operators are 
sometimes viewed with suspicion, as they are conceived 
to be fraught with uncontrolled approximations regard- 
ing the projection unto the physical sector of the Hilbert 
space necessary in those spin representations. Here we 
are using an exact method of projection onto the physi- 
cal part of Hilbert space that works even on the lattice 
(see below). 

The paper is organized as follows: Sec. |TT] introduces 
the model, the auxiliary-fermion representation and the 
projection schemes in detail. Simple Mean-Field approx- 
imations are discussed in Sec. IIIII where we demonstrate 
that these approaches are not able to capture frustra- 
tion effects but rather reproduce classical results. To this 
end in Sec. IIVI we introduce a phenomenological pseudo- 
particle lifetime that mimics quantum fluctuations. The 
results on the magnetization, susceptibilities and spatial 



spin correlations show that in a certain parameter range 
for this lifetime, the correct phase diagram is obtained. 
The main part of the paper, given by Sec. |V] is devoted 
to FRG. This method enables us to calculate the aux- 
iliary particle damping rather than treating it as an in- 
put of the approximation. In Sec. IV Al we first point 
out that the often applied static FRG scheme does not 
lead to a finite flow of the damping. Therefore in Sec. 
IV Bl in the framework of the standard truncation of the 
FRG equations we include the full dynamics. It turns 
out that within the latter (one-loop) approximation the 
strength of quantum fluctuations in the highly frustrated 
region is still underestimated, and a regime without Neel 
order or Collinear order is not found. We trace this de- 
ficiency of the one-loop approximation to the neglect of 
higher order contributions. Another way of saying this 
is that the Ward identities resulting from spin conser- 
vation are badly violated in the one-loop scheme, such 
that not even the RPA approximation is reproduced in 
that approximation. As shown by Katanin— the lat- 
ter problem may be remedied by using a dressed single 
scale propagator, thus including three-particle correla- 
tions with non-overlapping loops. As shown in Sec. IV CI 
using the Katanin truncation scheme we find a phase dia- 
gram in excellent agreement with results from numerical 
methods. In order to investigate the properties of the 
non-magnetic phase, correlation functions for columnar 
dimer and plaquette order are calculated in Sec. IV Dl 
We find that correlations for both kinds of dimerizations 
are clearly enhanced. Finally, the paper closes with a 
summary in Sec. I VII 

It is worth mentioning that although in this first 
presentation of our work using the newly developed 
FRG method we concentrated on demonstrating that the 
method is capable of giving results in agreement with 
results obtained mainly by purely numerical means, it 
should be clear that the method holds in fact considerable 
promise for future applications. First of all, it allows to 
treat thermodynamic, in contrast to finite size systems. 
Second, it is ready to calculate dynamical properties (at 
least on the imaginary frequency axis). Thirdly, it is eas- 
ily generalized to finite temperature (work in progress). 
Further, it allows in principle to address questions of crit- 
ical behavior near a quantum critical point. 



II. THE MODEL 

The effects of frustration in quantum spin models have 
been intensely studied in recent years. These models offer 
the possibility to investigate quantum phase transitions 41 
between magnetically ordered and disordered phases. Es- 
pecially in the context of deconfined criticality in two di- 
mensional spin system s 30 ' 42 , quantum phase transitions 
are the object of renewed interest. A standard model cap- 
turing these phenomena is the spin-I/2 Heisenberg model 
on a square lattice with an antiferromagnetic nearest 
neighbor coupling J\ > and a frustrating next-nearest 
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neighbor coupling J 2 > 0, see Refs. HEJIUllliilM 
H = J 1 J2s i -S j + J 2 J2^i-^j • (!) 

nn nnn 

As far as the ground state of the model is concerned, 
two limiting cases are well understood: For J 2 = the 
system is in a Neel-ordered phase with a magnetization 
of ~ 60% of the saturation magnetization. In the limit 
g = jj. — > oo the model reduces to two decoupled square 
lattices. Neel order on each of these lattices gives rise to 
the so-called Collinear long-range order with magnetic 
wave vectors Q = (ir, 0) or Q — (0, it). These two degen- 
erate ground states correspond to a parallel alignment of 
the spins along the y-axes and a antiparallel alignment of 
neighboring spins along the x-ax.es in the first case and 
vice versa in the second case. Increasing J 2 in the first 
limit or J\ in the second limit leads in both cases to frus- 
tration and to a decrease of the respective order param- 
eter, possibly all the way to zero. Indeed, the existence 
of a non-magnetic phase, indicated by numerical studies, 
approximately in the parameter region g « 0.4 . . . 0.65 
is widely accepte d 18 ' 20 ' 21 ' 28 ' 44 ' 46 . However, the nature 
of the magnetically disordered phase as well as the or- 
der of the quantum phase transitions is not known with 
certainty so far. Candidates for this phase are a spin 
liquid 7 » - 48 ' 49 and a valence-bond solid state (VBS). For 
the latter, different types of order have been proposed, 
e.g., columnar dime r 19 i 28 ' 46 i 51 and plaquett o 21 i 24 ' 28 i 47 or- 
der. Several studies give evidence that the transition 
from the non-magnetic phase to the Collinear phase is of 
first order— ~ 21 ' 28 ' 46 . We also mention that in the classical 
large spin-limit no magnetically disordered phase exists. 
Instead there is a direct first order transition between the 
Neel phase and the Collinear phase at g = 1/2. In this 
limit the respective magnetizations reach the saturation 
value. 

In the past the model has been studied with a 
variety of methods. Examples are: analytical ap- 
proaches based on field-theory method o 22 ' 29 or spin- 
wave theory 7 -^; numerical approaches such as exact 
diagonalizatio n 18 i 23 i 44 i 47 i 50 , coupled cluster metho d 20 ' 45 , 
series expansion method s 19 ' 28 ' 29 ' 46 ' 51 and Quantum 
Monte Carlo method 4 ^ . 

In this paper we address the ground state properties in 
a rather different way. In order to allow for application 
of Feynman-diagram technique s 52 ' 53 we reformulate the 
problem in a fermionic language by introducing auxiliary 
fermion o 31 ' 39 ' 40 . We represent the spin operators in terms 
of auxiliary fermions fi a , 

sf = ~E4<^ ■ (2) 

a/3 

Here cr' 1 (p, = x,y,z) are Pauli matrices, a, (3 =t,4 arc 
spin indices and i is the site index. We use units with 
ft, = fcs = 1. By construction, the representation ^ sat- 
isfies the correct commutation relations. However, the 
Hilbert space for a single site i is now spanned by four 



states, of which two, representing an empty and a dou- 
bly occupied site are unphysical. The projection to the 
physical sector of Hilbert space is given by the auxiliary- 
particle constraint 

Qi = £4j- = 1 (3) 

OL 

We present two different projection schemes to account 
for this constraint. 

A convenient approximate approach is to replace the 
constraint Qi — 1 by its thermodynamic average, (Qi) — 
1 . For a translationinvariant state the latter conditions 
are identical at each site, such that only a single condition 
remains. Since the constraint amounts to removing two 
of the four states per site, it is on average equivalent to 
half-filling of the system, which in case of particle-hole 
symmetry is effected by applying a chemical potential 
H = to the pseudofermion system. 

A different approach allowing for an exact treatment 
of the constraint even for lattice systems has been pro- 
posed by Popov and Fedotov 54 . It amounts to apply- 
ing a homogeneous, imaginary-valued chemical potential 
^ppv _ _ i^r_ ^ w j lere j 1 j s the temperature. Thus, within 
this scheme, the Hamiltonian H is replaced by 

H — > H ppv = H- /i ppv Qi ■ ( 4 ) 

i 

Note that H denotes the Hamiltonian (JXJ) using the repre- 
sentation of spin operators ((2|) . Given a physical operator 
O (i.e., an arbitrary sum or product of spin operators) it 
can be shown^ - that the expectation value (0) ppv , calcu- 
lated with H ppv and the entire Hilbert space, is identical 
to the physical expectation value (O) , where the average 
is performed with the original Hamiltonian H . The pro- 
jection works by virtue of a mutual cancellation of the un- 
physical contributions of the sectors Qi = and Qi = 2, 
at each site. It should be emphasized that although the 
Hamiltonian H ppv is no longer hermitian, the quantity 
(C) ppv comes out real- valued. If on the other hand O is 
unphysical in the sense that it is non-zero in the unphys- 
ical sector, e.g., the operator O = Qi, the expectation 
value (Qi) ppv is meaningless and one has (Qi) ^ (Qi) ppv - 
This approach is applicable to spin model o 40 ' 55 ' 56 but 
unfortunately it can not be extended to cases away from 
half filling. Although fj, ppv vanishes in the limit T — > 0, in 
principle the exact projection with /i — fi ppv and the av- 
erage projection with \i = are not equivalent at T = 0. 
This is due to the fact that the computation of an aver- 
age (. . . ) ppv does not necessarily commute with the limit 
T — > 0. Nevertheless it can be expected that in the model 
considered here both projection schemes are identical at 
T = 0. This can be understood with the following argu- 
ment: Starting from the physical ("true") ground state, 
a fluctuation of the pseudofermion number results in two 
sites with unphysical occupation numbers, one with no 
and one with two fermions. Since these sites carry spin 
zero the sector of the Hamiltonian with that occupation 
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is identical to the physical Hamiltonian where the two 
sites are effectively missing. Thus a fluctuation from the 
ground state into this sector costs the binding energy 
of the two sites which is of the order of the exchange 
coupling, even in the case of strong frustration- ^ 21 ! 24 ' 46 . 
Consequently, at T — pseudofermion-number fluctua- 
tions are not allowed and it is sufficient to use the simpler 
average projection with /i = 0. In most calculations we 
restrict ourself to this method. However, we again em- 
phasize that at T > both schemes differ. 

In the following we will formulate approximations in 
terms of resummed perturbation theory in the exchange 
couplings Ji , J 2 . The basic building blocks are the four- 
fermion interactions and the bare fermion Green's func- 
tion in real space 



1 



/< 



8ij6 a p ; fx 



inT 



(5) 



uj = (2n+ 1)ttT are the fermionic Matsubara- frequencies. 
Note that in diagrammatic expansions the Green's func- 
tions remain strictly local, i.e., Gy, a/ s = SijGi. a p. The 
momentum dependence in correlators like the suscepti- 
bility is generated by the non-local exchange couplings. 

We begin the calculations with a simple Mean-Field 
approach. It should be stressed that in our model a small 
parameter is absent. Accordingly, a controlled summa- 
tion of diagrams is a difficult task. For this reason we ex- 
tend the Mean-Field approach and set up a phenomeno- 
logical theory which explores the consequences of cer- 
tain assumptions on the width of the auxiliary fermion 
spectral-function and which gives qualitatively correct re- 
sults. 

Furthermore, the feasibility of diagrammatic approxi- 
mations allows the application of the Functional Renor- 
malization Group method (FRG)£Zr— . This scheme gen- 
erates an exact, infinite hierarchy of coupled differential 
equations for the one-particle irreducible m-particle ver- 
tex functions by introducing an infrared cutoff. In or- 
der to be able to solve these equations numerically one 
truncates the hierarchy of equations. The truncation is 
expected to give good results for not too strong interac- 
tion. It will turn out that the truncation procedure is a 
non-trivial problem for the model considered here. Ef- 
fectively FRG sums up infinite classes of diagrams. This 
is a crucial property in the present problem for which a 
small parameter does not exist. 

So far, FRG has been applied to low dimensional, inter- 
acting fermion systems, e.g., the two dimensional Hub- 
bard model§&£i, the single impurity Anderson model 62 
and the Luttinger liquid with impurities 6 ^, but a pure 
spin model has not been tackled with FRG. 



III. MEAN-FIELD THEORY 

The most elementary approximation for a spin model 
is the Mean-Field theory. In our fermionic description it 
corresponds to the Hartree approximation shown in Fig. 




FIG. 1: Diagrammatic representation of the Hartree approx- 
imation. The full line is the bare Green's function G°, Eq. 
([5]), the double-stroke line the self consistent one. The dashed 
line represents the interaction Ji or J2 and the dots are Pauli 
matrices x 1/2. 



[T] The closed loop of the renormalized propagator acts 
as the self-consistent Mean-Field. 

Note that the Fock term is exactly zero, since the non- 
local exchange coupling connects two points of the same 
fermion line. Dropping the requirement of exact pro- 
jection, one may allow for fermion hopping and make a 
Mean-Field ansatz with nonlocal propagators. The corre- 
sponding symmetry-broken phase is the so called resonat- 
ing valence bond (RVB)2 or the flux phas o 32 ' 64 . We will 
not consider Mean-Field amplitudes violating the auxil- 
iary particle constraint in this paper. 

By contrast, the magnetic order parameter (Si) = 

\ *52ap{fia a fip) that appears in the Hartree approxi- 
mation is a physical quantity. In the following calculation 
we also consider finite temperatures. Dyson's equation in 
Fig. [Dreads 



Gi{iuj) = [(iui + /Lt)l — £j(iw)]" 



(6) 



where G, S and 1 are matrices in spin space. The self- 
energy is coupled back to the renormalized Green's func- 
tion by 

£,M = \Y1 J ^ E a "\ E TrKG^o/)]^ . 



The couplings are written in the form , which is J\ if 
i,j are nearest neighbors, and J2 if i,j are next nearest 
neighbors. The factor e lw<5 , with an infinitesimal S > 0, is 
needed for the convergence of the Matsubara sum. If we 
assume magnetism along the z-direction, the self-energy 
has the form 



(8) 



To describe Neel- and Collinear order we split the lattice 
up into two sublattices A and B. In case of Neel order A 
and B form a staggered pattern while for Collinear order 
they form rows (or equivalently columns). Furthermore 
we require 



m = m ieA = -m ieB 



(9) 

Inserting Eq. © into Eq. (O and using Eq. ([SJ one 
obtains 



- a E J V a E E 



c 



8 ' iuj + ii — Crrij 



(10) 
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Using ± J2 



1 e*" 



/(z) and/(z-//PP v ) = SS5 l TT (/ is 



the Fermi function) one finds the following self-consistent 
equations for m for both types of order and both projec- 
tion schemes, 

M , , , /(Ji- J 2 )tanh(^) for/i = 

Neel-order: m = < , ' , ' , 

\ ( Ji - J 2 ) tanh(m/3) for /i = ^pp v 

(11a) 

_ ... , U tanh(^) for M = 

Collmcar-ordcr: m = < , ' 

1 J 2 tanh(m/3) for /z = ^pp v 

(lib) 

The spin polarization or, in short, magnetization Mj is 
given by 



M i = {St) = \^Tt[a*G i {iu)\e iuS 



(12) 



From the comparison of Eq. (O and Eq. (|12l) and using 
M^a = —Mi^B one finds a relation between and Mi, 



2M,(J 2 - 
-2M;J 2 



Ji) for Neel-order 

for Collinear order 
(13) 

From Eqs. (|lla[) and (11 lb)) the critical temperatures 
T c N6cl and T c Co1 can be determined. The instability with 
the larger transition temperature controls the type of or- 
der at a given g = This leads to 



< a < - : T r = T. N6cl = I 2 



4(l- 5 ) for M = 



Jl(l-fl) for M = M PP v 

(14a) 



a > I ■ T = T Co1 = J f ° r M = ° 

y ~ 2 ' c c \J ig for M = /iPP v 



(14b) 



Apparently, within this approximation, no non- 
magnetic phase is found at T = . Instead there is 
a first order transition from Neel- to Collinear order at 
g = i. The magnetization M = |M<|, which can be 
obtained from Eqs. (jllal) . (jl lb[) and (fl!?1) . reaches the 
saturation value M = \ at T = 0, and the classical large 
spin behavior is reproduced. These properties hold for 
both projection schemes. However, this is no longer the 
case for T > 0. The contribution of unphysical states 
with S = leads to a reduction of the magnetization 
in the average projection scheme. Also the critical tem- 
peratures come out a factor of two smaller in the aver- 
age projection scheme. The self-consistent equations for 
fi = jU ppv are identical to those obtained within the con- 
ventional Mean-Field theory in terms of spin operators, 
confirming that the cancellation of the unphysical states 
works correctly in this approximation. 

In summary, the simple Mean-Field theory leads to a 
Neel phase at g < \ and a Collinear-ordered phase g> \ 
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FIG. 2: Magnetizations MjMeei an d Mq q \ versus g within the 
phenomenological theory for different damping parameters 7 
(color online). 



but is not sufficient to describe the effect of frustration 
in destroying magnetic order in the regime g rss 

IV. FINITE PSEUDOFERMION LIFETIME 

In the Mean-Field approximation the effect of fermion 
scattering in generating a finite lifetime is not taken into 
account. We now briefly discuss a phenomenological 
framework for the ground state, which introduces the life- 
time r as a phenomenological parameter. To this end we 
model the retarded Green's function by, 

G r (uj) = — - — , E = -i-y with 7 = - . (15) 

CJ + 17 T 



The spectral function p(u>) acquires a finite width 7, 

= --Im G R (u) = -^7^ — 7 • (16) 



An analytic continuation of the self-energy E to the upper 
complex half-plane provides 



E(z) = —17 for Im z > 



(17) 



Since p(u>) is an even function it follows immediately 
from the spectral representation G{iuj) — de 
that G(z) and E(z) are odd functions with vanishing real 
parts along the complex Matsubara-axes. Thus we ob- 
tain 



G(iu) 



1 



iu> + 17 sgn(w) 



(18) 



To proceed, we need to specify the ^-dependence of the 
damping parameter 7. For J 2 = we put 7 in the form 
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FIG. 3: Phase diagram in the 7-g-plane. The dotted line 
shows the ^-dependence of 7 according to Eq. (1191) for 7 = 
0.36. 



7 = 7 J 1 where 7 is a dimensionless parameter. A similar 
situation is encountered for J\ — > 0, J 2 > 0, where the 
system is split up into two square lattices, each only with 
nearest- neighbor couplings Ji- Therefore in this limit 
the relation 7 = 7J2 holds. To interpolate between both 
limiting cases we assume 



7(^2) =lJi V1 + .9 2 



A. Hartree approximation 



(19) 



Replacing the bare Green's function by Eq. (fl8|) we 
now calculate the ground state magnetization within the 
Hartree approximation of Sec. IIIII In the limit T — > 0, 
using 4 J2iuj ~hi I ^ w ' tne new Mean- Field equation 
is given by 



4 ^ 



1 f°° 
/ dw £ 

Z7r J ~°° c=±i 



iu) + 17 sgn(w) — (rrij 



(20) 

Here it is obvious that the two projection schemes are 
identical because a shift of the Matsubara frequencies 
by fi ppv — — jjj becomes irrelevant in the limit T — > 0, 
provided that the Green's function, or equivalently the 
fermion spectral- function, is regular at u) = 0. Eq. ([20]) 
leads to the following self-consistent equations for the 
Neel- and Collinear magnetizations, 



M- 



M, 



1 /2M N6o l(Jl 

Neel = - arctan 

7T V 7 

1 /2McolJ2 

Col = — arctan 

7T V 7 



■/ 2 



(21a) 
(21b) 



The solutions of these equations are shown in Fig. [2] 
for different parameters 7. The case 7 = represents the 
Hartree approximation from Sec. IIIII An increase of 7 re- 
duces the magnetizations, especially in the region of high 
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FIG. 4: Self-consistent RPA equation for the susceptibility \ 
in diagrammatic representation. II denotes a single bubble. 



frustration. In particular, for small 7 there is still a direct 
transition between the two types of order at g = ^ , while 
for sufficiently large 7 a non-magnetic phase emerges. 
It appears that a broadening of the pseudofermion lev- 
els captures much of the effect of frustration expected to 
reduce or destroy magnetic order. In contrast to the sim- 
ple Mean-Field theory one now finds second order phase 
transitions and a Mean-Field critical exponent j3 = | of 
the magnetization. From the self-consistent equations a 
phase diagram in the 7-g-plane can be drawn, see Fig. [3] 
It shows only a narrow parameter range for 7 where the 
theory provides meaningful values for the phase bound- 
aries. For that reason it will be difficult to determine 
the damping parameter in approximative schemes. For 
example 7 = 0.36 leads to transitions at g c \ ks 0.39, 
g C 2 ~ 0.69 and also a realistic value for the magnetiza- 
tion at g = 0, i.e., MnooI ~ 0.35. This value of the width 
parameter 7 will be used in the following section to study 
the properties of the non-magnetic phase. 



B. Random Phase approximation 

In this section we calculate the spin susceptibility in 
the paramagnetic phase within RPA using the Green's 
function introduced in the last section. Fig. @] dis- 
plays the approximation in diagrammatic form. Since 
the RPA scheme can be obtained from the Hartree ap- 
proximation by taking the derivative with respect to 
the self-consistent field, the quantum phase transitions 
are located at the same point in both approaches. The 
conserving approximation scheme in the sense of Baym 
and KadanofT 65 "^ is an essential aspect here because spin 
conservation is an important constraint on the dynamics 
of the physical system; in addition the auxiliary particle 
constraint which allows only one particle per site requires 
a conserved particle number. In the non-magnetic phase, 
the equations in Fig. 0] are translation- invariant and can 
be transformed into momentum space. The static sus- 
ceptibility x(p, iv — 0) then has the form 



X(p,iv = 0) 



(U(w = 0))- 1 + J(p) 



(22) 



Here J(p) is the Fourier-transform of the interaction 
Jij = Jj_j which is given by 

J(p) = 2 Ji [cos(p x ) + cos(p y )] + 4 J 2 cos(p x ) cos(p y ) . 

(23) 
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FIG. 5: Static susceptibility for wave vectors (n, n) and (n, 0) 
and a damping parameter 7 = 0.36. The dashed lines visual- 
ize the phase boundaries (color online). 
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FIG. 6: Static correlation function x(R> w = 0) f° r distances 
R = |R| along a lattice direction, d is measured in units of a 
lattice spacing. Again 7 = 0.36 is used. In the upper panel g 
is slightly above the critical point g c \ and in the lower panel 
slightly below g C 2- 



Inserting the propagator from Eq. (fT5)) into a single bub- 
ble TL(iv = 0) this quantity is found as 



n(w = 0) = Ii zz {iv = 0) 
1 1 f . ( 1 



dw 



iuj + ij sgn(cj) 



Tr[(a 2 ) 2 



1 

27T7 

(24) 



The susceptibility in Eq. <J22]) together with Eqs. (J23J) , 
(fM| and (Unj) is evaluated for p = (n,ir) and p = 
(it, 0), the relevant wave vectors in the case of Neel- and 
Collinear order, respectively. The results are shown in 
Fig. [S] As expected for continuous phase transitions, the 
susceptibility with wave vector p = (it, n) (p = (0, 7r)) di- 
verges in the limit g c — > g c i + (g c g C 2 — 0). 



FIG. 7: Static correlation length £ in units of a lattice spacing 
versus g. The spectral width is 7 = 0.36. 



Finally, we discuss the static correlation function 
x(R,i^ = 0), which is obtained by transforming the sus- 
ceptibility from Eq. ([22|) into real space, 



X (R> - 0) = 



1 



(27T) 



dp x / dp. 



[n(^ = o)]-i + j( P )' 



(25) 

Evaluating Eq. (I2~5j) numerically with 7 = 0.36 for dis- 
tances R along the vertical or horizontal lattice direction, 
= i?e M , fj, = x,y, leads to the behavior shown in Fig. 
[6] For g slightly above the lower critical value g c \ (upper 
panel) the signature of the Neel phase is clearly seen. The 
correlation function forms a staggered pattern and the 
envelopes for positive and negative values only differ by a 
sign. At large enough distances R the data points are well 
fitted by an exponential decay while at small distances 
the decrease is faster. Inside the paramagnetic phase the 
envelopes are no longer symmetric around x(i?) = 0. For 
g slightly below the upper critical point g C 2 (lower panel) 
the correlation function still exhibits a staggered sign but 
the correlation between spins with an odd distance seems 
to vanish on approaching the critical point. Again, for 
large R an exponential function can be fitted and the 
correlation length is identical for even and odd distances. 
The asymmetry of the two envelopes can be understood 
by the fact that for Collinear fluctuations, two degen- 
erate patterns exist, the alignment of spins along rows 
and along columns. Thus near the upper critical point 
correlations are a superposition of both, 



X(R) = (-l) R die « +a 2 e 



(26) 



with (almost) identical weights a\ =0,2- Obviously this 
suppresses correlations for odd distances. Here £ denotes 
the correlation length. Away from the upper critical 
point Neel-like fluctuations emerge and we have a\ > 02- 
Eventually at the lower critical point a 2 vanishes. 

The correlation length £ is plotted in Fig. [7J The 
data indicate divergences at the phase boundaries but 
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get rather small in the vicinity of g — 0.4, i.e., down to 
£ pa 1.5. Remarkably, the smallest values for the correla- 
tion length are not reached at g = 0.5 where one would 
classically expect the strongest frustration. 

The phenomenological theory presented suggests that 
a broadening of the fermions' spectral function controls 
the phase diagram and the behavior of many physical 
quantities like the magnetization, the susceptibility and 
the spatial correlation function. However, a statement 
about the nature of the paramagnetic phase (columnar 
dimer or plaquette order) can not be made. Also critical 
behavior beyond Mean-Field is not accessible. Qualita- 
tively correct results are obtained by tuning the width 
7. However, this phenomenological parameter is not cal- 
culated within the theory and there is no simple way to 
calculate it. Unfortunately, summing up diagrammatic 
contributions of the Green's function to gain reasonable 
values for the width is a difficult tas k 39 ' 40 , because the re- 
sults strongly depend on the choice of the diagrams. For 
example, in the approximation of taking the self-energy 
to second order one gets a spectral function of the form 



(27) 



Then in the whole parameter range g > the gap A 
turns out to be so large that magnetic order is destroyed. 
On the other hand, from a completely self-consistent cal- 
culation of the second order self-energy, one finds that 
the spectral function is too narrow to allow for a para- 
magnetic phase. 

Note that there is no justification for a perturbative 
treatment in finite order. Instead diagram classes up to 
infinite order have to be summed. To approach this prob- 
lem in a more systematic way we will now use the Func- 
tional Renormalization Group method (FRG). 



V. FUNCTIONAL RENORMALIZATION 
GROUP METHOD 

The FRG method allows to sum up infinite classes of 
contributions in perturbation theory in a systematic way. 
So far, this method has been used to describe the weak 
coupling regime, e.g. of the Hubbard or Anderson mod- 
els. Our model does not have a small coupling constant, 
since there is no kinetic energy term in a spin Hamilto- 
nian. We nonetheless employ FRG in the usual way of 
neglecting higher order (three particle and higher order) 
correlation functions. As it turns out that this is not suf- 
ficient, we add higher order correlations in the form of 
self-energy corrections. Within FRG, one starts with the 
high energy sector, where Green's functions and coupling 
functions are known, and successively adds lower energy 
contributions. As a first step we define the cutoff proce- 
dure to be used by the following zero temperature bare 
Green's function, 

G 0A M = Q(H - A)G°M = 9(M ~ A) . (28) 

iu + [i 



d 

dA 1 



d 

dA 



2\W2' 2\ 



"X 1 h 



V' 



*V j V i fy A 



^4 3^ - ^3 4^ + ^4 3^ + ^ 3 4 ^ 
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FIG. 8: FRG equations for the self-energy and for the two- 
particle vertex. The line with an arrow is the full Green's 
function G A (iui) (see Eq. (|3ip ) and the line with an arrow 
and a slash is the single-scale propagator S A (iuj) (see Eq. 

El). 



In this cutoff dependent propagator all modes with \uj\ < 
A are projected out. For the rest of the paper we ap- 
ply the average auxiliary fermion projection scheme with 
\i = , as it is exact at zero temperature. However, 
it is not too difficult to implement the exact projec- 
tion scheme 5 ^, which increases the numerical effort by 
roughly a factor of eight. In the one- particle irreducible 
(1PI) version of FRG££~— employed here, G 0A (iuj) is in- 
serted into the generating functional of the 1PI vertex 
functions^in place of G°(iui). Taking the derivative with 
respect to A, an exact, infinite hierarchy of coupled dif- 
ferential equations for the vertex functions is obtained. 
To be more precise, the flow of the one-particle vertex, 
the self-energy E , depends on E and the two-particle 
vertex T. In turn, the flow of T depends on E, T and the 
three-particle vertex T^, and so on. At the end of the 
flow at A = when the theory is cutoff-free, the exact 
vertex functions are obtaine d 59 ' 60 . However, in explicit 
calculations one can only deal with a finite set of equa- 
tions and hence a truncation scheme has to be applied. 
Usually, by applying a weak coupling approximation, the 
three-particle vertex and higher vertices are neglected, 
resulting in a closed set of equations for E and T. This 
scheme will be applied in Sees. IV Al and IV Bl while in 
Sees. IV CI and IVDI we make use of an improved trunca- 
tion scheme^, which takes into account contributions of 
the three-particle type. For the conventional truncation 
scheme the equations for the self-energy E and the two- 
particle vertex T are depicted in Fig. [8]. In explicit form, 
these equations read 



A E A(i) = r A (l,2;l,2)S A M 



(29) 



dA 



r A (i', 2'; i, 2) = -L ]T [r A (F, 2'; 3, 4)r A (3, 4 ; 1, 2) 



3.4 

-r A (l', 4; 1, 3)r A (3, 2'; 4, 2)-(3o4) 
+r A (2', 4; 1, 3)r A (3, 1'; 4, 2)+(3o4)] 
G A MS A M (30) 
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Here the numbers are shorthand notations for the fre- 
quency, the site index and the spin index, that is 1 = 
{wi, ii, ai}, and J2i stands for an integral over u>i and 
sums over i\ and ct\. The full propagator G A (iu>) reads 



, A ,,, e(H-A) 



iuj — Y, A (iio) 

and the so-called single-scale propagator is defined by 



(31) 



^M = (g a m) 2 A (g oa m) - = _'5(I- 



A) 



(32) 

For the last expression in this equation a relationship 67 
for the product of 6-functions and (5-functions has been 
used. Note that G A and S A are local and translation in- 
variant in real space and proportional to the unit matrix 
in spin space. Thus the propagators in Fig. [S] and Eqs. 
([29]) and (|30|) carry only one composite index. 

Next we specify the initial conditions for the flow equa- 
tions at A = oo. In this limit, the free propagator van- 
ishes identically. Thus only one-particle potentials for 
the self-energy and bare interactions for the two-particle 
vertex remain. In the following we confine ourselves to 
the nonmagnetic phase. We defer consideration of the 
flow of the magnetic order parameter to later work. Ac- 
cordingly the free Green's function (|2"51) does not contain 
a one-particle field m that breaks rotational symmetry 
as in Eqs. ((6|) and (J5|). Note that although within this 
scheme the magnetic phases are not accessible, magnetic 
instabilities may be detected as divergences in the sus- 
ceptibilities. We have a vanishing self-energy for A = oo, 



^A— oo 



(iw) = 



(33) 



In this limit, the two-particle vertex is given by the bare 
interactions in antisymmetrized form, 



J- l 1 i' i 1 ^^ - J iiii 2"a,i(»i n Ojid! V'l 

— Jiii2 ~^ <J a 1 ,a 2 9 (T a 2 , cti <^V *2 ^i 3 > H 
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FIG. 9: Phase diagram in the 7-g-plane for a static FRG ap- 
proximation including a phenomenological parameter 7 (full 
line) . The dotted line shows the phase boundaries of the phe- 
nomenological theory from Fig. [3] 



particle vertex can now be represented as 

r A (i',2';i,2)= [ r^^^^^K^^ 
+ r<;i lll2 , W2' ; ^1, W 2 

+ Td hi 2 ( W l' I W 2' 5 U 2 , LO\ )S ai , a2 Sa 2 , ai ] 

^/ia^vu ■ (35) 

The indices s/d correspond to spin and density inter- 
actions. Note that energy conservation is implied, i.e., 
u)y + cl>2' = <jj\ + uj2- As another consequence of rota- 
tional invariance, the self-energy is an odd function in 
the frequency, as already pointed out after Eq. (TP?)) . In 
analogy to Eq. (fTT|) we write 



S A (iw) = -ij A (uj) 



(36) 



Inserting Eqs. (|3"T ]) . ([52 ]) . (|55 |) and ([55 ]l into Eqs. 
and ([3H)l the flow equations for 7, r s and Td can be cal- 
culated. 



(34) 



Here the factors \<J^a originate from the bare vertices 
and the Kronecker-<5 ensures that there is no fermion 
hopping on the lattice. Since the rotational invariance 
of the initial conditions is conserved during the flow, 
the two-particle vertex at hnite A is parametrized by 
spin-interaction terms oc c^cr^ and density-interaction 
terms oc S a pSjs- Since the propagators are local, the site 
index of an ingoing leg has to be identical to the site in- 
dex of the corresponding outgoing leg, which results in a 
total dependence on only two sites, i.e., i\ and To be 
more precise, translation invariance further reduces the 
site dependence only to the separation \i\ — Taking 
into account the antisymmetry in all variables the two- 



A. Static FRG 

Before considering the general case with all frequency 
dependencies, in this section we briefly discuss a static 
approximation 63 ' 68 . Putting all frequency arguments of 
the self-energy and vertex functions equal to zero leads, 
however, to a trivial solution since in that case the self- 
energy will be identically zero, provided it is assumed to 
be a continuous function of frequency. Therefore, in or- 
der to allow for a broadening of the spectrum we again 
assume the discontinuous form £ A = — i7 A sgn(w). How- 
ever, inserting this form together with the static two- 
particle vertex into the first flow equation (|2"9"|) leads to a 
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FIG. 10: FRG equation for the RPA scheme. The lines inside 
the boxes indicate the connections of the external legs. 



vanishing flow for 7 A due to the integration over an odd 
function on the right side. Obviously a static approx- 
imation can only be applied to the two-particle vertex, 
and 7 A has to be considered again as a phenomenologi- 
cal parameter that is independent of A. Using the static 
version of Eq. ([3"5]) . i.e., 

r (1)2; 1,2)= [r s i 1 i 2 Ga 1 ,a 1 cr a 2 ia2~^~^d ui 2 Q i ^ Q 2' a 2 ] 
Si 1 ,t 2 Si 2 , il (37) 

and the phenomcnological assumption E(ioi) = 
—ij sgn(oj) one obtains 



d 

dA' 



r A . . =- [Vr A r A - -2(r A 

sina 7r (A + 7) 2 l^-r 1 s 1 s J 12 s «i* 



s i±i2 \ s iiii d 



—r A = 2 1 [Vr A r A 

dA dllJ2 7r(A + 7 ) 2 L^ dllJ d " 

-r A f3F A -+- r A ") 

d iii 2 V s ~ d ) 



(38a) 



(38b) 



Note that the frequency dependence of E only affects the 
internal integration. By comparing Eq. (1341) and Eq. 
([57)) the initial conditions for T s and Yd can be read off, 



^A— oo 



4" 




7 Jii 



(39a) 
(39b) 



Solving Eq. (|38bj) with the initial condition (|39b[) gives 
r A - ia = 0. A finite set of equations for T A iii2 is obtained 
by neglecting all vertices with the distance \i± — ^| ex- 
ceeding a certain cutoff value. The resulting equations 
are solved numerically. A flow towards finite values for 
A — > indicates a paramagnetic phase while a diverging 
flow is a sign of a magnetic instability. The type of or- 
der can be extracted by transforming T A iii2 into Fourier- 
space and identifying the fastest momentum component. 
As a result one can draw a phase diagram in the j-g- 
plane, see Fig. [§] The figure compares the phase bound- 
aries with the results from the phenomenological theory 
in Sec. IIVI Again the boundaries are given by straight 
lines, but the value of the frustration parameter for which 
the paramagnetic phase has its largest extent moved from 



g = 0.5 to g ss 0.62. Interestingly, the phase diagram 
from Sec. IIVI can be reproduced within FRG. Consider 
the flow equation depicted in Fig. [TO] It only contains 
the second and third term of the right side of the equa- 
tion in Fig. [8] Furthermore it specifies how the ingoing 
and outgoing lines are connected. An examination of the 
terms in Fig. [8] reveals that the contributions in Fig. [10] 
are the only ones that couple two-particle vertices with 
different spatial separations of the outer legs among each 
other. This is due to the internal fermion bubbles in Fig. 
[TOl which can be located on an arbitrary site. For the 
other terms the vertex on the left side of the flow equa- 
tion is coupled only to itself or to the local vertex. The 
explicit equation corresponding to Fig. [10] reads 



— r 

dA 31112 



-,A 



1 



7T (A + 7) : 



S llj S J 12 



(40) 



again with a phenomenological 7. It can be show: 
that Eq. (j40|) reproduces the static RPA scheme and the 
phase diagram of Sec. IIVI Evidently the terms in Fig. 
[TOl are essential to obtain magnetism since they are the 
only ones that are able to describe collective phenomena. 
On the other hand, the remaining terms in Fig. [5] are 
only corrections that do not modify the phase diagram 
qualitatively. 



B. Dynamic FRG 



The considerations in Sec. IV Bl led to the conclusion 
that a static approximation of the FRG equations does 
not allow to calculate the central quantity governing the 
destruction of long range order: the pseudo fermion spec- 
tral width 7 . We will now treat the FRG equations in 
its full complexity and consider the dynamics with all 
frequency dependences as well as the back coupling of 
the self-energy into the two-particle vertex. This will 
lead to a finite spectral broadening without further as- 
sumptions. Again we make use of the truncation scheme 
that omits all vertices higher than the two-particle ver- 
tex. Inserting Eqs. j3T]), ([32]), ([35]) and (JSBJ into Eqs. 
([2"9"]l and (13T)]) . after a lengthy but straightforward cal- 
culation we end up with the flow equations and initial 
conditions presented in the Appendix. For convenience 
we write the two-particle vertex as a function of the in- 
variant frequency variables s, t, u, 



(41) 



defined by s = uii + LJ2, t = u)y — ui±, u — ujy — uji- The 
advantage of this parametrization is that r s and Td are 
both invariant under each of the transformations s — > — s, 
t —t, u — > —u, which can be deduced by a careful 
examination of the flow equations. This simplifies the 
numerics since only positive s, t, u have to be considered. 

In order to solve these equations numerically the con- 
tinuous frequencies will be discretized. We use a combi- 
nation of a linear and a logarithmic mesh. Since the two- 
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FIG. 11: Flow of the static Neel susceptibility for different 
frustrations (color online) 



FIG. 12: Flow of the static Collinear susceptibility for differ- 
ent frustrations (color online) 



particle vertex depends on three frequencies the compu- 
tational effort grows with the third power of the number 
of discrete frequencies. Regarding the truncation in real 
space the computing time grows with the fourth power of 
the length of the longest two-particle vertex (in two di- 
mensions, counting an internal site-summation, see Eqs. 
(TO) and dOb ). 

From the numerical solution we obtain physical quan- 
tities like the static correlation function Xiji^ = 0) by 
connecting the fermion lines of the two-particle vertex, 

/>oo 

= 0) = / dr(T T {Si(T)S*(0)}) 



iv=0 



A iv=0 



(42) 



To calculate this diagram two frequency integrals have 
to be performed numerically. The physical correlation 
function is recovered in the limit A = but we also con- 
sider Xij a t finite A. Transforming \% into momentum 
space, we obtain the magnetic susceptibility x A (p). The 
results are plotted in Fig. [TT] and Fig. [T2J It is seen 
that during the flow the Neel susceptibility (p = (jr, tt)) 
exhibits a divergence for all g < 0.55. On the other 
hand the Collinear susceptibility appears to diverge for 
all g > 0.55. In particular there is no parameter region 
without a magnetic instability and with a flow down to 
A = 0. In the present approximation the paramagnetic 
phase is obviously missing. The abrupt stop of the flow of 
the susceptibilities for g > 0.6 in Fig. [TTIand g < 0.55 in 
Fig. [12] can be traced to the divergence of the respective 
other susceptibility. For g at the transition, i.e., between 
0.55 and 0.6, the divergence is clearly indicated at the 
smallest accessible A. 

The absence of a nonmagnetic phase is quite unsat- 
isfactory. Obviously the spectral width comes out too 
small in this approximation. However, an essential im- 



provement is made in the next section where we use a 
different truncation scheme. 



C. Katanin truncation 

The only approximation that is involved in the FRG 
scheme described above is the truncation procedure in 
the hierarchy of differential equations. Unfortunately, 
the simple truncation employed above violates conserva- 
tion laws, expressed in terms of Ward identities. In order 
to improve the fulfillment of Ward identities Katanin de- 
veloped a one-particle self-consistent version of the two- 
loop FRG equations^. The basic modification there is 
the substitution of the single-scale propagator S A , see 
Eq. (|3"2"1) . by the total derivative of — G A with respect to 
A, 

S A (iuj) -> -^rG A M = 5 A M-(G A M) 2 ^-E A M . 
dA K J dA 

(43) 

It can be show n 37 ' 38 that such an approach is equivalent 
to an RPA + Hartree approximation if only terms of the 
RPA type (Fig. [TO]) are kept in the flow equation for 
the two-particle vertex. In this case Ward identities gen- 
erated by spin conservation are fulfilled exactly. As an 
application in a different context, for the reduced BCS 
model of superconductors exact Mean-Field results have 
been reproduced 38 . In particular, in conjunction with a 
small symmetry-breaking external field this scheme al- 
lows to access symmetry broken phase s 38 ' 69 ' 70 . If one 
keeps the terms additional to RPA on the right hand side 
of the second flow equation (see Fig. [FJ and also Ref. 
l6fj[ ). as we do, the exact conservation property is lost, 
but the remaining symmetry violating terms are gener- 
ated by overlapping loopdiagrams and may be expected 
to be smaller (see Ref. [33). While on the one hand the 
Katanin truncation scheme assures that the terms of an 
RPA+Hartree resummation are correctly included, we 
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find that the non-RPA terms are essential in providing 
just the right size of a finite auxiliary particle spectral 
line width. In that sense the non-RPA terms play a cru- 
cial role here: they control the pscudofcrmion damping 
and therefore the size of the nonmagnetic region in the 
phase diagram. 

As described in Ref. [H the substitution (|43| is made 
in Eq. but not in the equation for the self-energy, 

Eq. (|29|) . In the present work the above-mentioned small 
symmetry-breaking field is not applied. This would break 
the invariance of the two-particle vertex under s — s, 
t — > —t, u — > — u and would generate additional terms in 
the spin parametrization (|35[) . Effectively, with the sub- 
stitution (|4"5)l also contributions from the three-particle 
vertex are included. In the numerical implementation, in 
the equations for the two-particle vertex the internal bub- 
ble P£, n , Eq. (IA.5|) . is replaced by the modified bubble 
pA atj Eq. (|A.6|) . Due to the last term in Eq. (|A.6j) which 
does not contain a <5-function, the internal frequency in- 
tegration has to be performed explicitly. As a result now 
the computing time grows with the fourth power of the 
number of discrete frequencies. 

Typically we use 64 frequencies and discard all two- 
particle vertices with a spatial extent larger than 7 lat- 
tice spacings in each direction. Note that this trunca- 
tion corresponds to a system with 14x14 sites and pe- 
riodic boundary conditions because the longest bond in 
such a system extends over 7x7 sites. Exploiting lat- 
tice symmetries we end up with approximately 2.5 • 10 6 
coupled differential equations. The numerically deter- 
mined coupling functions and self-energies are inserted 
into Eq. (|42|) to calculate the susceptibilities shown in 
Fig. [T31 In the course of the flow the Neel susceptibil- 
ity for g = 0.2 shows a pronounced increase while the 
Collinear one stays very small, see Fig. [T^a). Obviously 
at that degree of frustration the system is in the Neel 
phase. However, we do not observe a real divergence of 
the susceptibility. When A gets too small the increase of 
the susceptibility stops and the flow exhibits an unstable 
and wiggly behavior which we attribute to numerical in- 
stabilities. Especially, for small A the flow is sensitive to 
the discretization of the frequencies. Going to larger sys- 
tem sizes the situation improves, i.e. one can follow the 
flow to larger susceptibilities and finds a steeper increase. 
Thus, unstable flows at small A can be identified as finite 
size effects. In the thermodynamic limit and with a suffi- 
cient number of discrete frequencies we expect a smooth, 
diverging solution indicating a magnetic instability. 

At g = 0.55, see Fig. [T5Tb). both susceptibilities ap- 
proach finite values for A — > 0, demonstrating the ex- 
istence of a phase with neither Neel nor Collinear long 
range order. Small oscillations are a consequence of the 
frequency mesh. In our numerics the limit A = can not 
be reached exactly because of an insufficient number of 
discrete frequencies at very low energy scales. Typically 
the flow is stopped at A « 0.0 Ui but can be easily ex- 
trapolated to A = 0. Finally at g — 0.8 the Collinear 
phase can be identified. In Fig. HUT c) the behavior 
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FIG. 13: Flow of the static susceptibility for the wave vectors 
p = (n,n) and p = (tt, 0) and different parameters g, (a) 
g = 0.2, (b) g = 0.55, (c) g = 0.8 (color online). 



is analogous to Fig. [T3Ta). but showing an increasing 
Collinear susceptibility. 

In order to investigate the properties of this phase fur- 
ther, we calculated the susceptibilities at additional pa- 
rameter values. The results in the physical limit A = 
are shown in Fig. 1141 Deep inside the paramagnetic 
phase our results are well-converged. With increasing g 
we observe a decreasing Neel susceptibility and an in- 
creasing Collinear susceptibility. The point where Neel- 
like fluctuations loose out compared to Collinear fluc- 
tuations lies at g « 0.6 in correspondence with the re- 



13 




(a) 



(b) 



■** «- 



FIG. 16: Patterns for (a) columnar dimerization and (b) pla- 
quette dimerization. Red bonds correspond to strengthened 
and green bonds to weakened interactions in the Hamiltonian 
H + F d or H + F p (color online). 



FIG. 14: Static susceptibility for the paramagnetic phase in 
the physical limit A = 0. The thick lines are obtained by a 
finite size scaling. The thin lines are the results for a maximal 
bond length of 7 x 7, 5 x 5 and 3x3 sites, from top to bottom 
(color online). 
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FIG. 15: Frequency dependent damping 7 A=0 (u;) obtained at 
the end of the FRG-flow compared to the constant damping 
with 7 = 0.36 used in Sec. IIVI (dashed line). Depicted is the 
case g — 0.55. 

suits in Sec. IV A\ i.e., clearly higher than the classical 
value g — 0.5. At the phase boundary to Collinear order, 
which turns out to be in the range g C 2 ~ 0.66... 0.68, 
the critical fluctuations require large system sizes, in or- 
der to obtain well-converged results. Here a finite size 
scaling (see thin red lines in Fig. [T4|) considerably en- 
hances the Collinear susceptibility and a beginning di- 
vergence is visible. The situation is very different near 
the phase boundary to Neel order. A divergence of the 
Neel susceptibility is not seen and finite size effects play 
a minor role. Instead, here the flow is highly sensitive 
to the frequency discretization, which causes large oscil- 
lations. Therefore it seems to be difficult to access this 
critical region and to obtain reliable data, see the dot- 
ted part of the blue line in Fig. [TJ] Here a denser mesh 



enhances and smooths the Neel susceptibility during the 
flow. An estimation of this phase boundary leads to the 
parameter region g c i w 0.4 . . . 0.45. 

Only few results on spin susceptibilities are found in 
the literature 1 ^. To the best of our knowledge so far 
no data are available for static spin susceptibilities. We 
note that our results on the phase boundaries are in good 
agreement with previous result o 18 i 20 ' 21 i 28 ' 44 i 46 . 

So far we notice that the behavior of the system near 
the two transitions is very different. To draw a conclusion 
concerning the order of the transitions, a closer investi- 
gation taking into account a flowing order parameter is 
necessary. 

Although the tendency towards formation of Neel- and 
Collinear phases has already been seen in the standard 
truncation of the previous subsection, the inclusion of 
certain higher order terms turns out to be essential. As 
shown above, the non-RPA like terms in conjunction with 
the Katanin truncation indeed lead to a damping 7 A (w) 
which is strong enough to generate a non-magnetic phase. 
The damping which is related to the self-energy via Eq. 
(|36| is no physical observable. This quantity is obtained 
from the first flow equation and can be compared with the 
frequency independent 7 used in the phenomenological 
theory from Sec. |lVl see Fig. [15] At g = 0.55, i.e., in 
the non-magnetic region we find that 7 A=0 (w) compares 
quite well with the choice 7 = 0.36 in Eq. (|T9| at relevant 
energy scales uj ~ J. In the region with high frustration, 
g 0.55, the static two-particle vertex T A iii2 (s = 0, t = 
0,u = 0), where i\, are nearest neighbors and using 
the initial conditions given by Eq. (IA.7I) . comes out as 
r A =° 2 (0,0,0) » 1.8 Ji at the end of the flow. 



D. Columnar dimer and plaquette order 

In this section we discuss the nature of the paramag- 
netic phase and investigate whether there is still some 
kind of long range order. Possible states currently under 
discussion are a spin liquid state (which does not break 
any symmetries) and a VBS state. For the VBS two 
dimerization patterns are of special interest: In a colum- 



14 




tion functions by 



FIG. 17: Main figure: Flowing dimer correlation Xa f° r 
g = 0.55. Inset: Dimer correlation Xd = ° m the limit A = 
versus g. 



nar dimer arrangement translation invariance along one 
lattice direction as well as rotation symmetry are bro- 
ken. For a plaquette valence-bond ordering, translation 
symmetry in both directions is broken while the rotation 
symmetry is intact. 

In order to probe the paramagnetic phase with re- 
spect to these states we add a small perturbation field 
to the Hamiltonian and investigate the response to 
1^19-21,23,28,29,47 In the context f FRG this concept has 

already been applied in Refs. [38.69. The fields can be 
chosen as 

F d = Sj2(-l)%, j S i+1 , j , 

f p = S J2 [(-^S^Si+U + (-l^'SijSij+j] , (44) 



Xd/p 



Jl 



((Si,j, Si+l,j)) d / p + ((Si+ij, Si+2,j)) d / p 



(45) 

Here the index d/p indicates that the correlator ((...)) is 



calculated with the Hamiltonian H + F, 



d/p- 



The factor 



S in the denominator eliminates the dependence on the 
strength of the perturbation such that x d / P start with 

1. An increase (decrease) dur- 



A=oo 



XT/ 



the initial value 
ing the flow shows that the system supports (rejects) the 
perturbation. Note that we again apply Katanin's trun- 
cation scheme. 

The columnar dimer correlation \d is plotted in Fig. 
[T7l It is seen that this quantity increases considerably 
during the flow. In the limit A = the perturbation 
is enhanced by a factor of w 3.8 but a divergence does 
not occur. We do not exclude an instability for this kind 
of order which might be masked here due to finite size 
effects. Remarkably we obtain plaquette correlations Xp 
with the same strength. Apparently the FRG scheme 
is not able to distinguish between dimer and plaquette 
correlations. This might be a consequence of the fact that 
the four-particle vertex from which such susceptibilities 
can be calculated directly, is not included in our FRG 
equations. The oscillations during the flow are again a 
consequence of the frequency mesh. 

Thus our results favor a spin liquid with enhanced 
equal time correlations Xd, Xp or a VBS with one of the 
two arrangements. Previous papers mainly calculated 
the columnar dimer and plaquette susceptibilities in the 
magnetically ordered phases rather than in the paramag- 
netic phasei^^^. 



VI. SUMMARY 



for the columnar dimer and plaquette order, respectively. 
Here i, j are components of the position vector and 6 is 
an energy much smaller than J\ and 3 2- Note that the 
expectation values (-Fd) and (F p ) are the ord er par ame- 
ters of these states. F d (see Refs. [TlllllIiii^E^ and 
F p (see Refs. 20.21) break the above-mentioned lattice 
symmetries and generate the two dimerization patterns 
shown in Fig. [TBI Possible instabilities should be visible 
as divergences in the corresponding equal time correla- 
tion functions Xd/ P = Jg^ 1 3=0- 

The coupling to these operators is included in the FRG 
formalism by modifying the initial conditions. The bare 
interactions in the limit A — > oo are slightly strengthened 
or weakened according to the dimerization patterns. Fur- 
thermore we have to take into account that due to the 
broken translation symmetries a two-particle vertex is 
no longer uniquely determined by one lattice vector. In 
the course of the flow we calculate the correlations of 
strengthened and weakened bonds and its relative differ- 
ence. We define equal time dimer and plaquette correla- 



The aim of this work is the development of new meth- 
ods and the calculation of properties for frustrated quan- 
tum spin models. We focused on the spin-1/2 J\-J% 
Heisenberg antiferromagnet on the square lattice, but 
our method is generally applicable to models of that 
type. Starting point of our approach is perturbation 
theory in the exchange couplings, summed to infinite 
order. In order to be able to use standard many-body 
techniques and to perform diagrammatic expansions, we 
applied the auxiliary-fermion representation of spin op- 
erators. To enforce the auxiliary-particle constraint, two 
different projection schemes have been employed: (1) en- 
forcement of the constraint on the average (which, how- 
ever, becomes exact at zero temperature), and (2) exact 
projection using an imaginary chemical potential^ . 

In a first exploratory study we use the RPA + Hartree 
approximation to access the ground state properties. 
Since the straightforward Mean-Field approach is not 
able to describe the suppression of magnetic order near 
g — y- — 0.5, we introduce a damping term in the bare 
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pseudofermion Green's function. This phenomenological 
self-energy accounts for scattering processes of the auxil- 
iary fermions which lead to a finite lifetime and a spectral 
broadening. We show that the damping reduces the mag- 
netic order especially in the regime of strong frustration. 
For sufficiently large damping one finds a paramagnetic 
phase around g = ^ = 0.5, as seen in numerical studies 
(we term this phase paramagnetic although it may pos- 
sess more complex magnetic correlations) . Furthermore, 
using RPA, we calculate the magnetic susceptibility and 
the spin correlations in the non-magnetic phase. We ob- 
serve critical behavior at the phase boundaries, i.e., a 
divergent susceptibility and correlation length. Within 
this method the basic properties come out qualitatively 
well, but a microscopic derivation of the pseudofermion 
damping is beyond the reach of simple diagrammatic re- 
summations. 

A more systematic approach is considered in the main 
part of the paper: the Functional Renormalization Group 
method (FRG). We employ its formulation in terms of 
the one-particle irreducible vertex functions and use a 
sharp frequency cutoff A. This method sums up large di- 
agram classes in a systematic way, e.g. the two-particle 
vertex function in the particle-particle channel and in 
two particle-hole channels and reaches therefore far be- 
yond the Mean-Field theory. Self-energy contributions 
are taken into account on equal footing with the vertex 
renormalizations. 

First we apply the conventional truncation scheme to 
the hierarchy of FRG differential equations, neglecting all 
vertices higher than the two-particle vertex. Further im- 
posing the static approximation, we find that the pseud- 
ofermion broadening is not generated in this way. Adding 
a phenomenological broadening, the results of the phe- 
nomenological theory are recovered, which is saying that 
RPA-likc diagrams in the particle-hole channel can be 
identified as the most important contributions. Includ- 
ing the frequency dependencies of the self-energy and the 
two-particle vertex we find magnetic instabilities in the 
whole parameter range. The latter is signalled by an 
immanent divergence of the susceptibilities at p = (tt, tt) 
and/or p = (tt, 0) in the course of the flow towards A = 0. 
It thus turns out that the truncation scheme is insuffi- 
cient to generate a strong pseudofermion damping and a 
nonmagnetic phase. This can be traced to the violation 
of Ward identities in this approximation. 

An improved approximation including self-energy ef- 
fects in the single-scale propagator has been suggested 
by Katanir>21. There the single-scale propagator is re- 
placed by the total derivative of the Green's function 
with respect to A. We find that the latter approach, 
even if only implemented on the one-loop level, repro- 
duces features of the Mean-Field theory and fulfills the 
Ward identities exactly, if the RPA-like contributions in 
the particle-hole channel are considered onl y 37 ' 38 . In our 
calculations we include the additional terms on the right 
side of the equation for the two-particle vertex in Fig. [5] 
Calculating the susceptibility, we are now able to distin- 



guish between the three phases. In particular we get a 
convergent flow down to A = in the region where the 
paramagnetic phase is expected. The phase boundary 
between Neel-order and paramagnetic phase is found to 
be at <7d « 0.4 . . . 0.45 and the transition from paramag- 
netism to Collinear order happens at g c i ~ 0.66 . . .0.68. 
Our findings agree well with the results obtained by other 
method a 18 ' 20 ' 21 ' 28 ' 44 ' 46 . Approaching the transition to 
Collinear order we observe a smooth divergence of the 
corresponding susceptibility. On the other hand near the 
transition to Neel-order a different picture emerges. Here 
it is difficult to access the critical region because the dis- 
cretization of the frequencies generates large oscillations 
in the flow at small A. 

Finally we probed the non-magnetic phase with respect 
to columnar dimerization and plaquette order by investi- 
gating the flow in the presence of appropriate small per- 
turbative fields. In the limit A = the correlations for 
both types of order are enhanced but a divergence is not 
found. These results indicate either strong dimer and 
plaquette fluctuations in a spin-liquid phase or a symme- 
try broken phase with dimer or plaquette order. 

The work reported here shows that in spite of the 
fact that quantum spin models are in the strong cou- 
pling regime by definition, partial resummations of per- 
turbation theory appear to capture the physics of frus- 
trated magnets at least on a qualitative level. The re- 
summations, done here in the framework of the Func- 
tional Renormalization Group method, account in a con- 
trolled and systematic way for all two-particle interaction 
processes, including all couplings between the different 
channels. If self-energy corrections are added in a consis- 
tent way, thereby including certain contributions of the 
three-particle vertex, it appears to be possible to calcu- 
late the central quantity of frustrated spin systems in the 
language of pseudoparticles, the damping of the pseud- 
ofermions in a controlled way. The systems we consider 
are in principle infinitely large but the spin correlations 
are only treated up to some finite length. Hence we do 
not have to deal with the effects of edges or periodic 
boundary conditions. However, the range over which our 
correlations extend typically include more than 200 sites 
which is much more than the system sizes accessible by 
exact diagonalization. Furthermore we do not make any 
assumption on the ground state or perform an expan- 
sion around any presumed state. Our starting point of 
free fermions without dispersion is completely featureless. 
Therefore our approach is straightforwardly applicable 
to a variety of models. Some obvious extensions of the 
present work are (1) consideration of other models such 
as the spin-1/2 J1-J2-J3 Heisenberg antiferromagnet on 
the square lattice, or geometrically frustrated models like 
the triangular or Kagome lattice, (2) generalization to 
finite temperature, (3) calculation of dynamical spin cor- 
relation functions. Work in this direction is in progress. 
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In this appendix we show the flow equations for 7, T s 
and Td with all frequency dependences, 



d A/ \ 1 

dK 1 M = Y 



- 2 (Id a (w + A, 0, w - A) - r£ tf ( w - A, 0, to + A)) 



+3(r* ii (w + A, w - A, 0) - r A 4 ( W - A, u + A, 0)) 
+r A „ (w + A, w - A, 0) - (w - A, w + A, 0) 



A + 7 A (A) ' 



(A.l) 



d_ 

dA 



= -f 



dw' 



(-2r^ li2 ( S; 



a/, wi' + w')r A - • (s, w 2 + a/, uii + oj') 



+ifi 1 i 2 (*»-w2' -w'.wi- W)r A fs^W,^ W) 



^A /„ , , , .' , . 1 , ,/\rA 



+r d 1 ili2 (s,-w 2 '-w',a; 1 ,+w')rsi 1 i 2 (s^2+w',a;i+a; / )) 
(P A (u/,s + o/) + P A (s + u;V)) 

-(2 ^ r A u ( Wl , + a/, t, wi - w')r* • l2 (W2 + w', t, -W2' + u/) 



S J l 2 

-A 



A 4ll2 (Wi/ + W', i, Wi - w'JI^fci! (W2 + W', -«? + W', *) 
*i lia K' + *> "I - "Old fci! ("2 + W', -W* + W', t) 

+r^ l4l + w', wi - w', i)r A iii2 (u 2 + w', t, -W2» + w') 
-r A uu (wi' + w', wi - w', t)r A iii2 (wa + w', i, -u 2 > + w')) 

(p A (u/,tw) + .p A (tw,u/)) 



S lil 2 \ 

+r A lll2 (^2' - w', -Wl - w', fi)r A 4l42 (W 2 - w', U>1' + w', u) 
+ r di 1 i 2 (w2' - u/,u)r A 4i42 (w 2 

(P A (u/, u + a/) + P A (w + a/, a/)) 



(A.2) 
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^ r dixi 2 (M> u ) 



(3r A li2 (s, -uj 2 ' - w'.wv + u')T^ lil2 (s, w 2 + w', u>! + J) 

+rdi li2 (s,-W2'-w',u;i/+w / )rd W2 (s^2+w' ) a;i+w')) 
(pV, s W) + P A (s W,w')) 

3 

-3I^ ilia (wi/ + a/, t, wi - w')r A Ml (wa + w', -W 2 , + u', t) 

-r* lll2 (Wl' + W', t, Wi - W')ld (^2 + W', -Wj- + W', t) 

-3r A Mi (w v + w', wi - w', t)r^ ixis (wa + w', t, -a* + w') 

("I' + W', Wi - W', t)Id *!* 2 (^2 + W', t, -W2/ + u/)) 

(P A (c/,fW) + -P A (*W,u/)) 
~( 3rA lll2 ( w 2' - w',u)r A 4i42 (w 2 -w'.wr 

+ r diii 2 ( w 2' -u/,-u>i -o/,u)r A ili2 (w 2 

(p A (o/, « + a/) + p a (u + ">'))' . 



(A.3) 



Note that the frequency parametrization of Eq. (HIT) is 
used for T s and r<j. The frequencies wi<, Wy, wi, w 2 on 
the right side stand for the inverse transformations 

Wl' = -{s + t + u) , W 2 ' = -{s-t-u) , 
wi = -(s-t + u) , w 2 = i(s + i-w) . (A.4) 

P A (a;i,w 2 ) denotes a bubble of 5 A and G A . For the 
conventional truncation as discussed in Sec. IVBI one gets 



- A) e(|wa|-A) 



Wi + 7 A (^i) w 2 + 7 A (w 2 ) 
(A.5) 

In this scheme the internal integration J do/ • • • simplifies 
to X^'=±a ' ' ' ■ For the Katanin truncation considered 



in Sec. IV CI we get a more complicated expression, 

5(|wi|-A) 9(|W2|-A) 



wi + 7 A (wi) w 2 + 7 A (w 2 ) 



_d A/ A e(K|-A) 9H-A) 

dA 7 1 lj y ( Wl + 7 A ( Wl )) 2 o; 2 +7 A ( W2 ) ' 1 ' 



In both cases P A (wi,o; 2 ) is an odd function separately 
in ui\ and w 2 . Finally from the comparison between Eq. 
([55)1 and Eq. ([55]) we get the following initial conditions, 

7 A=oo (cj)=0 , 

r A ,=£(M,*0 = Ri 2 , r A =? 2 (M,u) = o .(a.7) 
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